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Abstract 

In this case study, we model a planet’s magnetic and gravitational fields using 
spherical harmonic functions. As an exercise, we analyze data on the Earth’s 
magnetic field collected by NASA’s MAGSAT spacecraft, and use it to derive a 
simple magnetic field model based on these spherical harmonic functions. 

1 Introduction 

There are many times when it is useful to create a mathematical model of some 
physical phenomenon: that is, a set of mathematical equations that summarizes 
the results of many observations. For example, the Earth has a magnetic field 
similar to the magnetic field of a bar magnet. Suppose we wish to estimate the 
magnitude and direction of the Earth’s magnetic field at some specific location 
on the Earth’s surface. How would we do that? We could search past records 
for measurements made by various people, hoping to find some that are near the 
point of interest, then try to interpolate between the observation points. This 
method would be quite cumbersome, though — it would require sifting through 
thousands of observations made by many different observers at many different 
times, trying to find some appropriate observations from which to interpolate. 

A simpler method is to create a mathematical model that summarizes all 
the observations by fitting them to a set of mathematical equations. Once the 
model has been created, computing an estimate of the Earth’s magnetic field at 
a specific location is easy: just insert the latitude, longitude, and altitude of the 
location of interest into the equations, and out comes the magnetic field vector. 
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Such a mathematical model also makes it easy to look for trends in the data 
(the drift of the magnetic poles with time, for example). 

A similar method may be used for modeling the Earth’s gravitational field. 
By fitting many observations of the magnitude and direction of the gravitational 
acceleration to a set of mathematical functions, one may create a mathematical 
model of the Earth’s gravitational field that can be useful for applications such 
as high-precision orbit predictions. 


2 Spherical Harmonics 

For modeling the magnetic and gravitational fields of the Earth or other planets, 
it is customary to use special functions called spherical harmonics. In a sense 
these are two-dimensional counterparts of the sine and cosine functions used in 
Fourier analysis. Given a set of data defined over the surface of a sphere (such 
as the Earth), one can fit the data to a series of spherical harmonics in much 
the same way as one can fit data defined on a circle (i.e. periodic data with a 
period of 2i r radians) to a Fourier series. 

Spherical harmonic functions are actually complex- valued functions. Instead 
of using those directly, we’ll create a series using, separately, the real and imagi- 
nary components of the spherical harmonics, which are the two sets of functions 

cos (m</>) P ; m ( cos 9), 

sin (m</>) PJ 71 (cos 9) 

respectively. Here 9 and (f> are the usual polar and azimuthal angles (respec- 
tively) in spherical polar coordinates, l and m are integer indices (with m < l) 
of special functions Pf n { cos 9) called associated Legendre functions of the first 
kind. [1] The first few such Legendre functions (through l = 3) are shown in 
Table 1. 
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POINTER. Associated Legendre Functions. 

One of the tricky parts about working with spherical harmonic analysis is 
that there are a number of different normalization conventions for the associated 
Legendre functions in use, each of which gives rise to different leading coefficients 
for P™(cos0). The functions shown in Table 1 use the so-called Schmidt normal- 
ization convention, which is the one most commonly used in geomagnetism. [2] 
When working in other areas, you may encounter other conventions. The MAT- 
LAB function legendreO calculates associated Legendre functions, and includes 
an option that allows you to select among several different normalizations. 

Note also that the notation P ; m indicates a function with two integer indices, 
l and m; to is not an exponent. 


Table 1. Associated Legendre functions, P™( cos0). 

Pg (cos 9) = 1 
P° (cos 9) = cos 9 
Pf (cos 9) = sin 9 
P° (cos 9) = i (3 cos 2 6—1) 

Pi (cos 9) = V3 sin 9 cos 9 
P|(cos0) = i-\/3sin 2 0 
P° (cos 9) = | (5 cos 3 9 — 3 cos 9) 

Pi (cos 9) = \ \/6 sin 9 { 5 cos 2 9 — 1) 

Pi (cos 9) = \ a/ 15 sin 2 9 cos 9 
Pi (cos 9) = j y/lO sin 3 9 


3 Magnetic Field Models 

Now suppose we wish to model the Earth’s magnetic field using spherical har- 
monics. The Earth’s magnetic field is a vector field , meaning that there is a 
magnetic field vector associated with each point in space. That complicates our 
analysis a bit, since it would seem to mean that we have to fit separate series 
to each of the three components of the magnetic field. But there’s a simpler 
method: suppose we are only interested in modeling the magnetic field outside 
the Earth, due to electric currents inside the Earth’s core. Then we can assume 
there are no electric currents present in the region of space at which we are mod- 
eling the magnetic field. With these assumptions, it turns out that Maxwell’s 
equations of electromagnetism allow us to take the Earth’s magnetic field vector 
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POINTER. Spherical Polar Coordinates. 

You will sometimes encounter different conventions for spherical polar coor- 
dinates. As used here, 9 is the polar angle (measured from the +z axis), and (f> 
is the azimuthal angle (measured in the x-y plane counterclockwise from the +x 
axis). The (north) latitude is the complement of the polar angle, or 90° — 9 , while 
the (east) longitude is the same angle as <f>. 

The coordinate r is the radial distance from the center of the Earth. It is 
related to the altitude h above sea level by (approximately) r = a + h, where a is 
the mean radius of the Earth (6371.2 km). 


B to be the gradient of a magnetic scalar potential V : [3-5] 


B = -/i 0 W, 


(1) 


where the magnetic field B has units of teslas (T), the magnetic scalar potential 
P has units of amperes (A), and p,Q = 47r x 10~ 7 newtons per square ampere (N 
A” 2 ) is a constant called the permeability of free space. The symbol V is the 
gradient operator ; in spherical polar coordinates, VP is 


dv „ ldv „ 

VP (r, 0, 4>) = — e r + - — e# 
or r od 


1 dV . 

rsin# dcj) 


( 2 ) 


and e r , eg, and are unit vectors in the r, 9 , and </> directions, respectively. 
(See Pointer box on spherical polar coordinates.) Now all we need to do is fit the 
magnetic scalar potential P to a single spherical harmonic series, which looks 
like this: [2, 6] 


V (r, 9, </>) 



i 

(g™ cos mf) + h™ sin m0)P™ (cos 9) 

m=0 


(3) 


where a = 6371.2 km is the mean radius of the Earth, r is the radial distance 
from the center of the Earth (r > a since we’ve assumed we’re outside the Earth), 
and g™ and h™ are the expansion coefficients that we need to determine. The 
l = 1 terms represent the dipole component of the magnetic field, l = 2 the 
quadrupole component, l = 3 the octupole component, and so on. 


Activity 1. (a) Why are there are no l = 0 terms included in the summations 
in Eq. (3)? (b) Published tables of coefficients of g™ and h™ do not include any 
h® coefficients. Why would that be? 


In Eq. (3), the integer N is called the order of the spherical harmonic expan- 
sion, and determines how many terms there will be in the series. In general, the 
more terms, the more accurately the series will represent the data. However, at 
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POINTER. Units. 

When doing scientific calculations of this sort, it’s crucial that you pay proper 
attention to units of measurement. An easy way to avoid problems with units is 
to make it a rule that all variables in your computer program will always be stored 
in base SI units (kilograms, meters, seconds, amperes), and derived units based 
on them. For this project that means using these units for all calculations inside 
the program: 

• Magnetic field B in teslas (T). 

• Magnetic scalar potential V in amperes (A). 

• All lengths (r and a) in meters (m). 

• Magnetic dipole moment m in A m 2 . 

• All angles in radians (rad). 

• Coefficients g™ and h™ in teslas. 

If all inputs to your calculations are in base SI units, then the results of the calcu- 
ations will automatically be in base SI units also. Note, however, that when you 
look up coefficients g™ and h™ in the literature, they will be given in nanoteslas 
(nT). 

Remember to convert the MAGSAT magnetic field data to base SI units after 
reading them from the data file. 


some point the magnitude of the terms is about the same size as the measure- 
ment errors in the data, so it makes little sense to carry the series beyond that 
point. At the time of this writing, this series is typically carried out to N = 13 
for the full-scale geomagnetic field model. [7, 8] However, for the purposes of 
this case study, we’ll carry out the series to just N = 3 to make things simpler. 
This will still give a reasonably good model of the Earth’s magnetic field; it just 
won’t be as accurate as the TV = 13 model. 


4 Determination of Coefficients 

To determine the coefficients g™ and /ij", let’s apply the gradient operator 
(Eq. (2)) to the spherical harmonic series for the magnetic scalar potential V 
(Eq. (3)). Each of the components of the resulting vector will give one of 
the components of the Earth’s magnetic field vector. In geomagnetism, we 
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customarily define the three components of the geomagnetic field as: 


X = 

-B e = 

Mo(W) e 

(northward component) 

(4) 

Y = 

+7?0 = 

-Mo(VY )0 

(eastward component) 

(5) 

Z = 

— B r = 

Mo(Vh) r 

(downward component) 

(6) 


Let’s now calculate the partial derivatives of the potential V using Eqs. (2) and 
(3) and Table 1, calculating the series out to N = 2; this will give eight series 
terms for each component. The result is 




Y 


Z 


Mo- 


ldV_ 

de 


q 3 ^ Q' >J Q' >J 

77 sin Ogl + — cos <j> cos 9g\ H — -sin <f> cos 9h\ 

rytij rytj 


r 3 

3 CL^ q ,4 

2 ^4 sin ( 26, )^2 + ^3^4 cos (t> cos{29)g\ 


V3< 


+V / 3^- sin </> cos(29)h\ + — - — r cos(2</>) sin(20)g ; 


2 r 4 


+ 


sin(2</>) sin(26 ) )/i2i 


\/3 a 4 
2~ r 4 

1 dV 

^ 0 rsin9 d(f> 

q3 q3 qA 

Oa? ”1 — t sin <bg\ - cos <bh\ + OoS + V^— r sin 6 cos 9g\ 

1 y*o 4 * ^*4 " 

a 4 a 4 

— \/3— 7 cos <f> cos 0/i2 + %/3-j sin(20) sin 0g 2 


(7) 


-V3 — r cos(2</>) sin 0h 2 , 
r 4 

dV 

M0-5- 

or 

—2— cos 9g° — 2— cos </>sin 6*5^ — 2— sin 0 sin 


-^-x( 3cos26 '- 1 )52- 


3V3a 4 


2 r 4 
3^ a 4 


3^ a 4 


cos ^>sin(20)(72 


sin <^sin(20)/i2 — 
sin(2</>) sin 2 0h 2 . 


3 v^a 4 


cos(2 (j>) sin 2 0^2 


( 8 ) 


(9) 


Activity 2. Continue calculating derivatives to expand each of these series to 
order N = 3. This will mean seven additional terms (g°, g 3, h\, <?f, h\, g|, 
and /i|) for a total of 15 terms for each of the X, Y, and Z components. (You 
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may wish to check your results using a symbolic mathematics program such as 
Mathematical 

As you can see, the number of terms increases rapidly as the order N is 
increased. For a spherical harmonic expansion of order N, each of the X , Y, 
and Z series will have (N + l) 2 — 1 terms. For the full N = 13 model, that’s 
195 terms in each series. That’s why we’re limiting the expansion to IV = 3, for 
which there are just 15 terms in each series. 


Observations of the Earth’s magnetic field will take the form of a set of com- 
ponents (X, Y, and Z), measured at a specific polar angle 9 (the complement of 
the latitude), longitude </>, and distance from the center of the Earth r. Our data 
set will then be a set of these variables, and we wish to solve for the unknown 
g™ and /i™ coefficients. 

One fairly straightforward way to do this is to write Eqs. (7) - (9) in matrix 
form. Doing this to order N = 2, we have Eq. (10): 
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Here there are M different data points, each of which consists of an observed 
magnetic field vector (Xj, k), and Z{), and the location r,;, 9i , and <pi at which 
that observation was made. The left-hand side is a 3 M x 1 column vector, 
and contains the observed magnetic field components; let’s call this vector b. 
The matrix on the right of the right-hand side is an 8 x 1 column vector that 
contains the coefficients we want to solve for; let’s call that vector g. The large 
matrix on the left of the right-hand side is a 3 M x 8 matrix that contains each 
of the terms of the spherical harmonic expansions of X , Y, and Z , evaluated 
at each location; we’ll call this matrix A. The first row of A times vector g 
gives the terms of the spherical harmonic expansion for X, evaluated at the 
location of data point 1; the second row times g gives the terms for Y evaluated 
at the location of data point 1; and the third row times g gives the terms for Z 
evaluated at data point 1. Rows 4-6 repeat the pattern: they’re the terms in 
the expansion for X, Y, and Z, but evaluated at the location of data point 2; 
rows 7-9 are evaluated at the location of data point 3, etc. Then we can write 
Eq. (10) concisely as 

b = Ag. (11) 

We know all the magnetic field observations (vector b), and all the spherical 
harmonic terms evaluated at the data locations (matrix A); we just need to 
solve for the coefficient vector g. Formally , we could do this with a matrix 
inverse, just by multiplying each side on the left by A . 

g = A : b. 

But there’s a problem here: matrix A isn’t necessarily square. In fact, it will typ- 
ically have far more rows than columns, and thus constitute an overdetermined 
system (i.e. there are more equations than unknowns). How can we compute the 
inverse of a non-square matrix? Technically we can’t, but there is a technique 
available that allows us to solve Eq. (11) for vector g in a least-squares sense: 
it returns the g that best fits the data in matrix A. This technique is called 
singular value decomposition. 


5 Singular Value Decomposition 

The singular value decomposition (SVD) of a matrix A expresses A as a product 
of three matrices: 

A = U£V t , (12) 

where in our case A is the matrix of spherical harmonic expansion terms (which 
is known), and all the matrices on the right-hand side are determined by the 
SVD method. Singular value decomposition exposes the geometric structure of 
a matrix, which is an important aspect of many matrix calculations. A matrix 
can be described as a transformation from one vector space (e.g. vector Ag) 
to another (e.g. vector b); the components of the SVD quantify the resulting 
change between the underlying geometry of those vector spaces. 
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(13) 


Substituting Eq. (12) into Eq. (11), we have 

b = (USV T )g. 

In Eqs. (12) and (13), matrices U and V are square and orthogonal, in the 
sense that their column vectors are orthonormal (i.e. U T U = I and V T V = 
I, and therefore U” 1 = U T and V” 1 = V T ). The matrix E is a generally 
non-square matrix containing what are called the singular values of matrix A 
(denoted cq) along its main diagonal and zeroes elsewhere. Once the singular 
value decomposition is performed (i.e. matrices U, X, and V have been found), 
then we can solve Eq. (13) for g: 

g = VX+U T b. (14) 

Here X + = (X T X) -1 X T is called the pseudo-inverse of matrix X [9,10]. 


Activity 3. Derive Eq. (14) from Eq. (13). 


Because of the structure of X (a non-square diagonal matrix), X + can be 
found by simply transposing X and taking the reciprocals of the elements along 
its main diagonal, so X+ = 1/cq, X+- = 0 for * ^ j, and X + X = 1. If X is 
non-singular (i.e. none of the eq are near zero), you don’t need to do any more. 
However, if X is singular or nearly so (i.e. if any of the cq are near zero) you will 
want to replace the small values of eq (those below a specified tolerance) along 
the main diagonal of X by setting the corresponding diagonal elements of X + to 
zero; this prevents small measurement errors in b from producing large changes 
in g. The replacement of the small values in X by zeros in X + is equivalent 
to introducing a smallest cutoff singular value cq along the main diagonal of 
X in order for the matrix to be invertible. However, the choice of such cutoff 
singular values is not unique. As a working solution, the optimal choice of such 
small singular values can be easily implemented by allowing a threshold ratio 
between the smallest and the largest singular values which is not smaller than 10 
times the machine precision. If this threshold is not met, you increase the cutoff 
singular values by 10 until it meets the criterion. In all cases you should be able 
to test the stability of the pseudo-inverse by changing this cutoff singular value. 

Eq. (14) is then the solution to our problem. Given the magnetic field 
observations in vector b and the singular value decomposition of matrix A 
(which gives U, V, and X), we have the desired coefficients g. But how do 
we perform the singular value decomposition of A? The simplest and most 
reliable method is to use published algorithms [11-15] or the function svd() in 
MATLAB; the details of the internal workings of these implementations are a 
bit too complex to go into here. 


Activity 4. Equation (14) is most efficiently calculated by doing the multiplica- 
tions from right to left. Verify this for yourself by calculating the total number 
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of floating-point multiplications required to compute this equation when the 
right-hand side is evaluated from right to left, and when it is evaluated from 
left to right. 

Hint: In Equation (14), if M is the total number of data points, and the 
model is of order N = 3, then what are the sizes of matrices V, E + , U T , and 
of vector b? 

Hint: When multiplying an to x n matrix by an n x p matrix, the total 
number of floating-point multiplications required is to x n x p. 


Challenge 1 . Implement Eq. (14) to order N = 3 using MATLAB or some 
other programming language to find the coefficients <7™ and h,™ of the Earth’s 
magnetic field. For data, you can use a set of observations made by the 
MAGSAT spacecraft during the years 1979-80. MAGSAT collected a magnetic 
field data point every 0.5 second for 6 months, for a total of about 30 million 
data points. To make the data a bit more manageable, we’ve selected roughly 
1 point out of every 300,000 (about 1 point every 40 hours) for a total of 100 
data points randomly distributed around the globe, and made this available as 
file magsat-data . dat. The MAGSAT data is in the form of a plain ASCII text 
file, in the following format: 


Column 

Description 

1 

Time of observation (milliseconds of day; ignore) 

2 

Latitude 90° — 6 (north positive; deg) 

3 

Longitude <j> (east positive; deg) 

4 

Radial distance from center of Earth r (km) 

5 

Magnetic field component X (nT) 

6 

Magnetic field component Y (nT) 

7 

Magnetic field component Z (nT) 

8 

Attitude processing flag (ignore) 


Your program will need to read this data, do any needed unit conversions, 
compute and store the matrix elements, and then solve the least-squares problem 
to find the coefficients g™ and h™. 

If you are using MATLAB, you may wish to use the built-in function svd() 
to perform the singular value decomposition. You may also wish to investigate 
using the “economy” version of the function, which takes the form svd(A,0) 
and eliminates unneeded rows of zeros to save space and computation time. 

Once you have your program working, you may compare your results with 
the first 15 coefficients of the DGRF 1980 magnetic field model, available from 
NO A A at http : / /www. ngdc . noaa.gov/IAGA/vmod/igrf . html. Finding the er- 
ror (percent difference) between each of your coefficients and the corresponding 
DGRF 1980 coefficient will give you a rough idea of how well you were able to 
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Figure 1: Map of the magnitude of the geomagnetic field at the surface of 
the Earth for 1980, produced from a spherical harmonic model. Contours are 
labeled in units of nanoteslas (nT). Note that the field is generally strongest 
(red) at the magnetic poles and weakest (blue) near the equator. Note also the 
large dip in magnitude around central South America and the South Atlantic 
Ocean; this is known as the South Atlantic Anomaly. 


model the geomagnetic field using just 100 points of Magsat data. Of course, 
your results will not match DGRF 1980 exactly, because you’re not using the 
same set of observations as input — you’re only using 100 points of MAGSAT 
data. Nevertheless, the comparison will be close enough for you to tell whether 
you’ve done the calculation correctly. 


Challenge 2. As an additional project, you may wish to try implementing 
a geomagnetic field model. This would mean writing a computer program to 
calculate the X , Y, and Z components of the Earth’s magnetic field vector at 
any given latitude, longitude, and radial distance from the center of the Earth, 
using Eqs. (1) through (3), along with the g f 1 and hf 1 coefficients you just 
calculated in the previous Challenge. If you have the program calculate the 
magnetic field vector at many points along a grid over the surface of the Earth 
(say, at 1° intervals in latitude and longitude, and assuming a spherical Earth 
of constant radius r = a = 6371.2 km), you can produce contour maps of each 
component. Figure 1 shows something similar: a contour map of the magnitude 
of the magnetic field vector (i.e. the square root of the sum of the squares of 
the A, Y, and Z components) at the Earth’s surface. 
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POINTER. The Inverse Tangent Function. 

The inverse tangent function in Eq. (18) may look a bit odd. Why not just 
cancel the two minus signs? The reason is that the signs of both the numerator 
and the denominator must be used to ensure that the longitude will be in the 
correct quadrant. 

When computing the inverse tangent of a ratio like this, the rule is: do the 
division, then compute the arctangent of the quotient. The result, as computed by 
a calculator or computer, will be between —90° and +90°. If the denominator of 
the original ratio was positive, then use this returned answer; if the denominator 
was negative, then add 180° (n radians) to the returned answer. In Eq. (18) we 
have to keep the minus signs as they are to ensure that this extra 180° is added 
when the denominator (including the minus sign) is negative. 

Most computer programming languages include a special built-in function to 
handle this case, usually called something like atan2(y,a;). This will compute 
the arctangent of y/x, then check the sign of x to place the result in the correct 
quadrant. Be sure you remember this when computing the arctangent in Eq. (18). 


6 Analysis 


Besides using them in the magnetic field model, the coefficients g™ and h r l n may 
also be used to directly calculate a few quantities of interest. For example, the 
magnetic dipole moment of the Earth’s magnetic field can be shown to be given 
by [16] 


m = 


4na 3 

Mo 


V / M ) 2 + (si 1 ) 2 + (M) 2 - 


(15) 


It turns out that the magnetic poles of the Earth are not located at the 
rotation poles; they are rather located some distance away, and move from one 
year to the next. We can compute some quantities related to the location of 
the magnetic poles directly from the coefficients g] 71 and /i™. For example, the 
tilt angle a of the magnetic axis relative to the rotation axis can be shown to 
be [17] 


cos a = 


vW+W+TMF 


(16) 


The geographic latitude ipjv and longitude A at of the magnetic north pole are 
found from [1] 


-ST 

sin (Djsr = — . : 

V(9°i) 2 + (g\) 2 + (h\) 2 

(17) 

, -hi 

tan Ajv = r 

~9 1 

(18) 
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Activity 5. Based on your earlier results where you found the magnetic field 
coefficients g™ and h™, determine (a) the magnetic dipole moment of the Earth; 
(6) the tilt of the magnetic axis; and (c) the latitude and longitude of the 
magnetic north pole. (All values you calculate will be valid as of 1980, when 
the MAGSAT observations were made.) 


7 Conclusions 

In this case study, you have learned how to model the Earth’s magnetic field us- 
ing MAGSAT data, spherical harmonic functions, and computational techniques 
for solving linear systems. Such models enable you to estimate the Earth’s ge- 
omagnetic field vector at any given location without having measurements of 
these vectors at that particular location. 

Using the techniques described here, you may also model a planet’s gravita- 
tional field. [18, 19] In that case, the measurements are of the acceleration due 
to gravity g, which can be written as the gradient of a gravitational potential 

g = - VS (19) 

The gravitational potential Q is expanded in a spherical harmonic series, just 
as was done with the magnetic scalar potential V. 

Spherical harmonic expansions have also found applications in plasma physics. 
It has recently been shown [20] that if you fit a plasma’s distribution function 
to a spherical harmonic series, then the moments of the plasma (i.e. the plasma 
density, bulk velocity, temperature, and pressure) can be quickly computed as 
functions of the expansion coefficients. This technique is similar to formulas 
for calculating the Earth’s magnetic dipole moment and magnetic pole location 
described earlier. 


About the Authors 

David G. Simpson is a physicist at the NASA Goddard Space Flight 
Center, where he has worked on flight software for the Hubble Space Tele- 
scope, as well as data analysis for instruments on the IMAGE and Cassini 
missions. He is also an adjunct professor of physics at Prince George’s Com- 
munity College in Largo, Maryland. He has a Ph.D. in applied physics from 
the University of Maryland at Baltimore County. He may be contacted at 
David.G.Simpson@nasa.gov, or visit http://caps.gsfc.nasa.gov/simpson. 

Adolfo F. Vinas is a senior staff research physicist at the NASA Goddard 
Space Flight Center, where he has worked on data analysis, theory and simu- 
lation of solar and astrophysical plasmas. His major research interests include 
studying plasma kinetic and magnetohydrodynamic (MHD) physics of the solar 


14 


wind, heliosphere, shock waves, MHD and kinetic simulation of plasma insta- 
bilities and turbulence associated with space, solar and astrophysical plasmas. 
He has a Ph.D. in physics from the Massachussetts Institute of Technology. 


References 

[1] J.R. Wertz. Spacecraft Attitude Determination and Control. D. Reidel, Dor- 
drecht, Holland, 1985. 

[2] J.A. Jacobs. Geomagnetism, Vol. 1. Academic Press, Orlando, Fla., 1987. 

[3] J.R. Reitz, F.J. Milford, and R.W. Christy. Foundations of Electromagnetic 
Theory, Third Edition. Addison- Wesley, Reading, Mass., 1979. 

[4] Officer, Charles B. Introduction to Theoretical Geophysics. Springer- Verlag, 
New York, 1974. 

[5] Garland, George D. Introduction to Geophysics. W.B. Saunders, Philadel- 
phia, 1971. 

[6] S. Chapman and J. Bartels, Geomagnetism (2 vol.). Oxford, London, 1940. 

[7] IAGA Working Group V-MOD. International Geomagnetic Reference Field: 
the eleventh generation. Geophys. J. Int. (2010) 183 , 1216-1230. 

[8] International Geomagnetic Reference Field, at http://www.ngdc.noaa. 
gov/I AGA/ vmod/ igrf . html 

[9] Moore, E.H. On the reciprocal of the general algebraic matrix. Bull. Amer. 
Math. Soc., 26 , 9, 394-395 (1920). 

[10] Penrose, Roger. A generalized inverse for matrices. Proc. Cambridge Phil. 
Soc., 51 , 406-413 (1955). 

[11] G. Strang, Introduction to Linear Algebra, Fourth Edition. Wellesley Cam- 
bridge, Wellesley Mass., 2009. 

[12] G.E. Forsythe, M.A. Malcolm, and C.B. Moler, Computer Methods for 
Mathematical Computations. Prentice-Hall, Englewood Cliffs N.J., 1977. 

[13] G.H. Golub and C.F. Van Loan, Matrix Computations, Third Edition. 
Johns Hopkins, Baltimore, 1996. 

[14] LAPACK — Linear Algebra PACKage, available at http://www.netlib. 
org/lapack/ 

[15] ARPACK software routines to solve large-scale eigenvalue problems, avail- 
able at http://www.caam.rice.edu/software/ARPACK/ 

[16] R.T. Merrill, M.W. McElhinny, and P.L. McFadden, The Magnetic Field 
of the Earth. Academic Press, New York, 1996. 


15 


[17] M.G. Kivelson and C.T. Russell. Introduction to Space Physics. Cambridge, 
New York, 1995. 

[18] M. Bursa and K. Pec. Gravity Field and Dynamics of the Earth. Springer- 
Verlag, Berlin, 1988. 

[19] W. Torge. Geodesy, Third Edition, de Gruyter, Berlin, 2001. 

[20] A.F. Vinas and C. Gurgiolo. Spherical harmonic analysis of particle veloc- 
ity distribution function: Comparison of moments and anisotropies using 
Cluster data. J. Geophys. Res., 114 , A01105, doi:10.1029/2008JA013633 
(2009). 


16 


